/*
 * UserDefinedFunc.cpp
 *
 *  Created on: Nov 11, 2011
 *      Author: yi
 */

#include "mbuild/DefaultFunctions.h"


Gaussian::Gaussian(void) {
	/*
	 * Construct the Guassian distribtuion function.
	 *
	 */
	  IsAdditive=true;
	  offset=0;
	  solute_value=Inf;
	  density=0;
	  std=1.0;
	  s2=std*std;
	  s1=density*0.39894228/s2;
	  mode=SOLUTE_NEQL;
}

Gaussian::~Gaussian() {

}

void Gaussian::set_width(double w)
{
	/*
	 * Gaussian distribution width. Equal to standard deviation of the sample.
	 */
	std=w;
	s2=2*std*std;
	s1=density*0.39894228/std;
}



void Gaussian::set_offset(double o) {
	/*
	 * Set the offset of the plane to the surface.
	 */
	offset=o;
}

void Gaussian::set_density(double density)
{
	this->density=density;
	s1=density/sqrt(2*PI*std*std);
}

void Gaussian::set_sigma(double sig)
{
	this->density=sig;
	s1=density/sqrt(2*PI*std*std);
//	printf("Sigma: %e\n",s1);
}

double Gaussian::f(double d)
{
	/*
	 * Return the Gaussion distribution.
	 */
	d-=offset;
//	printf("s1:%f",s1);
	return s1*exp(-d*d/s2); //s1=density/sqrt(2*PI*std*std);s2=2*std^2;
}

void Gaussian::Solute_Do_Not_Eql(double v) {
	/*
	 * Solute should not have value equal to v.
	 */
	mode=SOLUTE_NEQL;
	solute_value=v;
}

double Gaussian::get_density()
{
	return density;
}

bool Gaussian::IfEqualSolute(double v)
{
	if(v != solute_value) return true;
	return false;
}

Switch::Switch() {
	/*
	 * Switch is additive
	 */
	IsAdditive=false;
	solute_value=Inf;
	in=-Inf;
	out=Inf;
	on=0;
}
Switch::~Switch() {

}
void Switch::set_values(double in, double on, double out)
{
	this->in=in;
	this->on=on;
	this->out=out;
}


void Switch::set_density(double d)
{
	//do nothing
}

double Switch::get_density()
{
	return 0.0;
}

void Switch::set_width(double w) {
	width=w;
}



double Switch::f(double d) {
//	if(d<0) printf("%f \n",d);
	if ( d>width) return (out);
	else if (d>0) return (on);
	else return(in);
}


Gouy::Gouy()
{
	mode=DONOTHING;
	  solute_value=Inf;
	  IsAdditive=true;
	  offset=0.0;
	  sigma=0.0;
	  rho=0.1;
	  z=1.0;
	  temp=298.0;
	  setup();
}

Gouy::~Gouy()
{
}



void Gouy::set_lipid(double area, double anfr,int charge)
{
	sigma=charge*anfr*Electron/area;
	setup();
}



void Gouy::set_ion(double valence, double conc)
{
	rho=conc;
	z=valence;
	setup();
}



void Gouy::set_temp(double temp)
{
	this->temp=temp;
	setup();
}





void Gouy::set_offset(double o)
{
	offset=o;
}

void Gouy::setup(void)
{
	phi0=-2334.2*sigma/sqrt(temp*rho);
	k=5.622667*z*sqrt(rho/temp);
	alpha=log(phi0+sqrt(phi0*phi0+1))/z;
	phi0=1.7235e-4*temp/z*log(phi0+sqrt(phi0*phi0+1));
	alpha=(exp(alpha)-1)/(exp(alpha)+1);
}
double Gouy::get_density()
{
	return 0.0;
}

void Gouy::set_density(double density){
	this->sigma=density;
	setup();
}


double Gouy::f(double d)
{
	double phi;
	d-=offset;
	if (d<=0)
	{
		return phi0;
	}
	else {
		phi=alpha*exp(-k*d);
		phi=(1.7233e-4)*temp*log((1+phi)/(1-phi));
	}
	return phi;
}

Nsmooth::Nsmooth() {
	IsAdditive=false;
	solute_value=Inf;
	nsmth=0;
	thick=0;
	in=0;
	out=1;
}



Nsmooth::~Nsmooth()
{
}



inline double Nsmooth::f(double d)
{
	double dt=d*2/thick+1;
	double t=1;
	if(dt > 99) return out;
	if(dt < 1e-2) return in;
	for(int i=0;i<nsmth;i++) {t*=dt;}
//	printf("%f %f \n" ,d,(in+(out-in)*t/(1+t)));
	return (in+(out-in)*t/(1+t));

}



void Nsmooth::set_density(double density)
{
	return;
}

void Nsmooth::set_density(double in,double out)
{
	  this->in=in;
	  this->out=out;
}

void Nsmooth::set_thick(double thick)
{
	if(thick<=0) throw ParameterError();
	this->thick=thick;
}



double Nsmooth::get_density()
{
	return 0.0;
}

void Nsmooth::set_nsmth(int n)
{
	this->nsmth=n;
}



Sigmoid::Sigmoid(){
	  IsAdditive=false;
	  solute_value=Inf;
}



Sigmoid::~Sigmoid(){
}



inline double Sigmoid::f(double d)
{
	  return out+(in-out)/(1+exp((d-alpha)/beta));
}



void Sigmoid::set_width(double w)
{
	  this->beta=w/4.394;
}



void Sigmoid::set_offset(double o)
{
	  this->alpha=o;
}
void Sigmoid::set_density(double d){
	return;
}


void Sigmoid::set_boundary(double b1, double b2, double width)
{
	 this->beta=width/fabs(log((in-b1)/(b1-out))-log((in-b2)/(b2-out)));
	 this->alpha=beta*log((b1-out)/(in-b1));
}



double Sigmoid::get_density()
{
	return 0.0;
}

void Sigmoid::set_density(double inside,double outside)
{
	  this->in=inside;
	  this->out=outside;
}
CylindricalHole::CylindricalHole()
{
	setup(0,0,0,0);
}

CylindricalHole::CylindricalHole(double r_up, double r_bottom, double length,double origin)
{
	setup(r_up,r_bottom,length,origin);
}



void CylindricalHole::setup(double r_up, double r_bottom, double length,double origin)
{
	if(r_up< 0 || r_bottom <0 || length <0) throw ParameterError();
	this->r_up=r_up;
	this->r_bottom=r_bottom;
	this->length=length;
	this->origin=origin;

}


double CylindricalHole::f(double x, double y, double z)
{
	if(length==0 ) throw NotInitialized();
	double r=(x*x+y*y),r0;
	r0= r_bottom+(z-origin)*(r_up-r_bottom)/length;
	if (z<origin || z> origin+length) return 1;
	if (r>r0*r0) return 1;
	else return -1;
}

CylindricalHole* CylindricalHole::clone()
{
	CylindricalHole* h=new CylindricalHole;
	h->setup(r_up,r_up,length,origin);
	return h;
}


















